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Abstract — The alternating direction method of multipliers 
(ADMM) has sparked recent interest as an efficient optimization 
tool for solving imaging inverse problems, such as deconvolution 
and reconstruction. ADMM-based approaches achieve state-of- 
the-art speed, by adopting a divide and conquer strategy that splits 
a hard problem into simpler, efficiently solvable sub-problems 
{e.g., using fast Fourier or wavelet transforms, or proximity op- 
erators with low computational cost). In deconvolution problems, 
one of these sub-problems involves a matrix inversion {i.e., solving 
a linear system), which can be performed efficiently (in the 
discrete Fourier domain) if the observation operator is circulant, 
that is, under periodic boundary conditions. This paper proposes 
an ADMM approach for image deconvolution in the more 
realistic scenario of unknown boundary conditions. To estimate 
the image and its unknown boundary, we model the observation 
operator as a composition of a cyclic convolution with a spatial 
mask that excludes those pixels where the cyclic convolution 
is invalid, i.e., the unknown boundary. The proposed method 
can also handle, at no additional cost, problems that combine 
inpating (recovery of missing pixels) and deblurring. We show 
that the resulting algorithm inherits the convergence guarantees 
of ADMM and illustrate its state-of-the-art performance on non- 
cyclic deblurring (with and without inpainting of interior pixels) 
under total-variation (TV) regularization. 

Index Terms — Image deconvolution, image deblurring, alter- 
nating direction method of multipliers (ADMM), variable split- 
ting, non-cyclic deconvolution, inpainting, total variation. 



I. Introduction 

The alternating direction method of multipliers (ADMM), 
originally proposed in the 1970's ifTSll . emerged recendy as a 
highly flexible and efficient tool for addressing several imaging 
inverse problems, such as denoising ll22ll . llT4l . deblurring 
12 1, inpainting |2|, reconstruction |18|, motion segmentation 
||26I, to mention only a few classical problems. ADMM-based 
approaches make use of variable splitting, thus allowing a 
straightforward treatment of various prior/regularization terms 
OJ, such as those based on frame representations and on the 
well known total-variation (TV) norm [[20l|, as well as the 
seamless inclusion of several types of constraints, such as 
positivity [22 J . 11 4i . ADMM is equivalent or closely related 
to other techniques, namely the so-called Bregman and split 
Bregman methods lfT6l . ||6l, ll25l and the Douglas-Rachford 
splitting Is), in, ini, mi. For a recent review of ADMM and 
its many appUcations, as well as a comprehensive literature 
review, see fS]. 

Most ADMM-based algorithms for imaging inverse prob- 
lems involve, as one of their steps, the solution of a linear 
system (equivalently, a large matrix inversion) 16), lT4ll . 
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112311 . This fact is simultaneously a blessing and a curse: on 
the one hand, the matrix being inverted is related to the 
Hessian of the underlying objective function, thus providing 
the algorithms with second order information and, arguably, 
justifying the excellent speed of these methods; on the other 
hand, this inversion (due to its typically huge size) limits their 
applicability to problems where it can be efficiently computed. 
For ADMM-based image deconvolution algorithms f6l, 
fT4|. l23j . this inversion can be efficiently computed using 
the fast Fourier transform (FFT), as long as the convolution 
operator is cyclic/periodic (or assumed to be so), thus diagonal 
in the discrete Fourier domain. However, in most real prob- 
lems, the convolution is not cyclic, thus the required matrix 
inversion can not be solved efficiently using the FFT. 

In an image deconvolution problem, the pixels located near 
the boundary of the observed image depend on pixels that 
lie outside of its domain; accordingly, it is necessary to make 
some assumptions about the values of those unobserved pixels, 
i.e., to adopt some boundary condition. The usual approach is 
to assume one of the classical boundary conditions imported 
from the partial differential equations literature, namely: peri- 
odic boundary condition; Dirichlet boundary condition, where 
the external boundary pixels are zero-valued; Neumann (or 
reflexive) boundary condition, where the external boundary is 
a reflection of interior region next to it; anti-reflexive boundary 
condition, where the external boundary is the negative of what 
it is in the Neumann case |8|, As illustrated in Figure [T] 
all these standard boundary conditions are quite artificial and 
do not correspond to realistic imaging systems; they are merely 
motivated by computational convenience considerations. 

Assuming a periodic boundary condition has the advantage 
of allowing a very fast implementation of the convolution as 
a point-wise multiplication in the discrete Fourier transform 
(DFT) domain, with both the direct and inverse DFT being 
very efficiently implemented by the FFT. However, in real 
problems, there is no reason for the external (unobserved) 
pixels to follow periodic (or any other) boundary conditions. 
A well known consequence of this mismatch is a degradation 
of the deconvolved images, namely the appearance of ringing 
artifacts emanating from the boundaries (see Figs. |2] (a) and 
[3] (a)). These artifacts can be reduced by pre-processing the 
image to reduce the spurious discontinuities at the image 
boundaries created by the periodicity assumption; this is 
what is done by the well-known "edgetaper" function in the 
MATLAB Image Processing Toolbox. In a more sophisticated 
version of this idea that was recently proposed, the observed 
image is extrapolated to create a larger image with smooth 
boundary conditions ITtII . 

Convolution under Dirichlet boundary conditions can still 
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Fig. 1. Illustration of the assumptions underlying the periodic, Neumann, 
and Dirichlet boundary conditions. 

be carried out in the DFT domain, using the classical zero- 
padding technique and embedding the non-periodic convolu- 
tion into a larger periodic one. However, in addition to the 
mismatch problem pointed out in the previous paragraph (there 
is usually no reason to assume that the boundary pixels are 
zero), this choice has another drawback: in the context of 
ADMM-based methods, the zero-padding technique does not 
allow the required matrix inversion (mentioned above) to be 
efficiently computed by resorting to the FFT |i4j. 

A. Contributions 

Instead of choosing a specific boundary condition or ap- 
plying a boundary smoothing strategy, a more realistic op- 
tion is to consider that the external boundary pixels are 
unknown, which is what actually happens in most imaging 
systems. This approach can be formulated as a simultaneous 
deconvolution and inpainting problem, where the unobserved 
boundary pixels around the image are estimated while the 
image is simultaneously being deconvolved. Although this 
formulation has been used before (see 1 8 1, ||T9l ), it has not been 
addressed using state-of-the-art ADMM-based methods. In 
fact, ADMM is particularly suited to handle this formulation. 



The observation model involves the composition of a periodic 
convolution with a spatial mask; the former corresponds to a 
point-wise multiplication in the spatial domain, while the later 
corresponds to a point-wise multiplication in the DFT domain. 
We will show how, in the context of the ADMM algorithm, 
these two operations can be decoupled and handled separately. 

Very recently, a related approach (although not expUcitly 
described as such) was proposed in flT], by exploiting ear- 
lier work for deconvolution without boundary artifacts under 
quadratic regularization lfT9ll . That approach requires solving, 
at each iteration, a linear system of the size of the number 
of unknown pixels, which can only be done numerically. 
In contrast, every update equation of the ADMM algorithm 
herein proposed has closed form, and the method is shown to 
satisfy sufficient conditions for convergence. 

The proposed method will be experimentally illustrated 
using total-variation (TV) deblurring (by adapting the state- 
of-the-art algorithm described in |23|); the results show the 
advantage of our approach over the use of the "edgetaper" 
function (in terms of improvement in SNR) and over an 
adaptation of the recent strategy of 121] (in terms of speed). 
Finally, our approach is also tested on inverse problems where 
the observation model consists of a (non-cyclic) deblurring 
followed by a generic loss of pixels (inpainting problems). 
Arguably due to its complexity, this problems has not been 
often addressed in the literature, although it is a relevant one: 
consider the problem of deblurring an image in which some 
pixel values are unaccessible, e.g., because they are saturated, 
thus should not be taken into account. 

B. Outline 

The remaining sections of the paper are organized as 
follows. Section HIl reviews the ADMM. Section Hill describes 
our ADMM-based approach in the context of TV-based image 
deconvolution. An experimental comparison of our method 
versus the use of the well known "edgetapper" function and 
the use of the recent strategy of I.21J , is reported in Section 
II VI that section also illustrates the use of the proposed method 
for simultaneous deblurring and inpaiting. Finally, Section |V] 
concludes the manuscript. 

Notation: we use bold lower case to denote vectors (e.g. 
X, y, u), and bold upper case (Roman or Greek) to denote 
matrices (e.g., A, D, T, F). The superscript ^ denotes 
vector/matrix transpose, or conjugate transpose in the case of 
complex vectors/matrices. 

II. The Alternating Direction Method of 
Multipliers (ADMM) 

A. The Standard ADMM Algorithm 

Consider an unconstrained optimization problem 

min/(z)+g(Gz), (1) 

where /:]&"■-;>» = »□ {+cx.} and g : Rp -> t ai-e proper, 
closed, convex functions, and G S R^^" is a matrix. The 
ADMM algorithm for this problem (presented in Algorithm[T]l 
can be motivated and derived from different perspectives (e.g.. 
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augmented Lagrangian, Douglas-Rachford splitting); for a 
recent comprehensive review, see |5 |. Convergence of ADMM 
was shown in |12j, where the following theorem was proved. 

Theorem 1: Consider problem ([T]l, where G e M^^" has 
full column rank and / : R'' ^ R and 5 : Rp ^ R are 
closed, proper, and convex functions. Consider arbitrary ^ > 
0, uo, do e W. Let r]k > and pk > 0, for fc = 0, 1, be 
two sequences such that J^'kLo Vk < 'X) and J2T=o < 
Consider three sequences S R",Ufe g R^*, and e MP, 
for fc = 0, 1, satisfying 



Zfc+i - argmin/(z) 



Ml 



Gz-Ufe-dfe||^ 



Ufe+i - argmin5(u) + -||Gzfe+i-u 
u 2 



< 



< 



Pk 



dfc+i dfc - (Gzfc+i - Ufc+i). 

Then, if ([T]i has a solution, say z*, the sequence {zfc} con- 
verges to z* . If ([T]i does not have a solution, then at least one 
of the sequences (ui,U2, ...) or (di,di, ...) diverges. 

Algorithm 1: ADMM 

1 Initialization: set fc = 0, choose /i > 0, Uq, and dg. 

2 repeat 

Zfc+i ^ argmin^/(z) + §|| G z - Ufc - dk\\l 
Ufe+i ^ argmin„.g(u) + f || G Zk+i - u - dfe||| 



^ dfe - (Gzfe+i - Uk+i ) 



fc fc + 1 
7 until stopping criterion is satisfied 



Under the conditions of Theorem 1, convergence holds for 
any value of /i > 0; however, the choice of this parameter may 
strongly affect the convergence speed |5|. It is also possible 
to replace the scalar p hy a positive diagonal matrix T = 
diag(/ii, pp), which means replacing the quadratic terms of 
the form ^||Gz-u-d||2 by (Gz - u - d)^T(Gz - u- d). 

B. The ADMM For Two or More Functions 

Following the formulation proposed in O, lfT4l . consider a 
generalization of ([T]i, with J > 2 functions, 

.7 

mm^5,(HWz), (2) 



where gj 



H 



IPjXn 



— R are proper, closed, convex functions, and 
are arbitrary matrices. We map this problem 



into the form ([T]i as follows: let /(z) = 0; define matrix G as 



where p = pi 



- pj, and let /2 : 



(3) 



R be defined as 
5(u)=5].g,(uW), (4) 



where each uW) G RPj is a sub-vector of u = 
[(uWr,...,(u(^)r]'^. 

The definitions in the previous paragraph lead to an ADMM 
for problem (|2]i with the following two key features. 

> Firstly, the fact that /(z) = turns line 3 of Algorithm[T] 
into an unconstrained quadratic problem. Letting T be a 
positive block diagonal matrix (which associates possibly 
different parameters pj to each function gj) 

PJ elements 

(5) 

the solution of this quadratic problem is given by (denot- 
ing C = Ufc + dk) 

argmm(Gz-C)^T(Gz-C) = (g^Tg) ^G'^TC 
J2p, (H(^))^Hwj"'x: A'^- (H(^'^)V^^ (6) 



T = diag( pi,...,pi , 

" V ' 

pi elements 



, Pj, pj^ , 
Pj elements 



A3) 



where C , u 



(j) „(i) 



and d 



with C'^' - "fe -1 , S , "fe , -fc 

for j = 1,..., J, are the sub-vectors of ^, u^, and d^, 
respectively, corresponding to the partition in (O, and 
the second equality results from the block structure of 
matrices G (in (O) and T (in (|5]l). 
Secondly, the separable structure of g (as defined in 
(01) allows decoupling the minimization in line 4 of 
Algorithm [T] into J independent minimizations, each of 
the form 



"i+i ^ arg min^ ^^(v) 



Pj_ 
2 



,0) 



|2 
I 2' 



(7) 



for j = 1,..., J, where s^^' = H^-J^Zfe+i - d)^'. This 
minimization corresponds to the so-called Moreau prox- 
imity operator of gj/ p (denoted as prox^.y^.) (see ^9), 
iTOl and references therein) applied to s^^^, thus 



arg mm — x- 
X 2 " 



=0) I 



p, 



For some functions, the Moreau proximity operators are 
known in closed form |9|; a well-known case is the 
£1 norm (7||x||i — 7X]i for which the proximity 
operator is the component-wise application of the soft 
threshold function: soft(w, 7) — sign(w) max{|ti| — 7, 0}. 
The resulting ADMM for problem (|2|l is presented in 
Algorithm |2l The convergence conditions given in Theorem 1 
for the standard ADMM translate, in the present formulation, 
to the following: (i) for 5 to be a proper, closed, convex 
function, it suffices that each function gj be itself proper, 
closed, and convex; (ii) for matrix G, as given in (|3]l, to 
be of full-column rank, a sufficient condition is that one the 
matrices H^^^ has full column rank. If G has full column 
rank, a sufficient condition for the inverse in ^ to exist is 
that Pj > 0, for j = 1, J. 

If the proximity operators of the functions gj are compu- 
tationally inexpensive, the computational bottleneck of Algo- 
rithm I2] is the matrix inversion in line 6. As shown in 1231 
and ||24) . these matrix inversions can be efficiently tackled 
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Algorithm 2: ADMM-2 



1 Initialization: set k — 0, choose ^i, > 0, uq, do- 



2 repeat 



C ^ Ufc + dfc 



Zfe+l 



(i) 



for j 1 to J do 

ui^,^prox^^./,^(H(^)z,+i-dif')) 

di!),^di^)-(H(^)z,+i-ui!^,) 

end 

k + 1 

10 until stopping criterion is satisfied 



using the FFT, if all the matrices H*^^^ are circulant, thus 
diagonalized by the DFlQ (more details will be given in the 
next section). Because of this condition, the work in f23\ 
and Il24 i was limited to deconvolution with periodic boundary 
conditions, under TV regularization. More recent work in I2j|, 
ID, IH has shown how this inversion can still be efficiently 
handled using the FFT (and other fast transforms) in other 
problems, such as image reconstruction from partial Fourier 
observations and inpainting, and with other regularizers, such 
as those based on tight frames. In the next section, we will 
show how to handle deconvolution problems with unknown 
boundary conditions. 

III. Proposed Approach 

This section presents the proposed approach for image de- 
convolution with unknown boundary conditions. The approach 
will be described in the context of image deconvolution under 
TV regularization ll20ll . ||23]| . as this is, arguably, the most 
standard regularizer for this class of imaging inverse problems. 
Extension to other regularizers, namely frame-based analysis 
or synthesis formulations fS], f3l, fT?|, is straightforward. Be- 
fore addressing the non-cyclic deblurring case, we review the 
cyclic deblurring case, which will provide the basic elements 
and structure for the proposed method. 



where: x,y e M" are vectors containing all the pixels 
(in lexicographic order) of the sharp and degraded images, 
respectively; A e R"^" is a matrix representation of the 
convolution (under periodic boundary conditions) with some 
kernel; e M'^^" is the matrix that computes the horizontal 
and vertical differences at pixel i (also with periodic boundary 
conditions) ||231 . Il24l ; A > is the regularization parameter. 
The sum V]"-] ||Dix||2 = TV(x) defines the so-called total- 
variation^ (TV) of image x. The first term in (O results from 
assuming that y is a noisy version of the convolved image 
Ax, where the noise is white and Gaussian with unit variance 
(assuming unit variance implies no loss of generality, since any 
other value may be absorbed by the regularization parameter 
A). 

Clearly, problem ^ has the form (O, with the following 
mapping: J = n + 1, 



9t 



gt(v) = A ||v||2, for i = 1, 



-i(v) = ^||y- v| 



gn+i ■■ K'" ^ K, 

h("+i) = a e m"^". 



(9) 
(10) 

(11) 
(12) 



Concerning matrix Y (see (|5]l), we adopt fij ~ a > 0, for 
j = 1, ...,n, and fi„+i = /3 > 0. With these choices, and the 
mapping in (l9Tl-(fT2]l. the fundamental steps of Algorithm |2] 
are as follows. 



As shown in 11231 . 1124 1. the matrix to be inverted in line 
4 of Algorithm |2] can be written as 

n 

a^DjDj +/3A^A = 

a(D'')'^D'' + a(D'')^D^ +/3A^A, (13) 

where D'' e M"^" (respectively, D" e M"^") is the 
matrix collecting the first row of each of the n matrices 
Di (respectively, the second row of each of the n matrices 
Di); that is, D'' computes all the horizontal differences 
and D" all the vertical differences. Assuming periodic 
boundary conditions, all the matrices in (fT3T l are circulant, 
thus can be factorized as A = U^AU, D'* = U"^ A'' U, 
and D" — A" U, where U is the unitary matrix 
(U"^ = U^) corresponding to the DFT, while A, A'\ 
and A" are diagonal matrices. Consequently, the inverse 
of the matrix in (T3[ can be written as 



A. TV-Based Deconvolution with Periodic Boundary Condi- 
tions 

The classical formulation of TV-based image deconvolution 
consists of a convex optimization problem 



arg mm 



1, 



xGR" 2 



Axil 



Aj^IlD, 



X 2, 



(8) 



1=1 



'in fact, since we are dealing witli two-dimensional (2D) images, these 
matrices need to be block-circulant with circulant blocks (BCCB) to be 
diagonalized by the 2D DFT. For simplicity, we refer to BCCB matrices 
simply as circulant and to the 2D DFT simply as DFT. 



U^(a|A'*P + a|A"|2 + /3|A|2)"'u, (14) 

which involves a diagonal inversion (with 0(n) cost) 
and the products by U and (the direct and inverse 
discrete Fourier transforms), which have 0{nlogn) cost, 
if implemented by the FFT. 
• The sum defining the vector to which this inverse is ap- 

-In fact, this is the so-called isotropic discrete approximation of the 
continuous total-variation norm Q, 1201 . 
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plied (line 4 of Algorithm |2]i can be expressed compactly 
as follows: 



J 

E 



(n+l) 



T/.(n+l) 



a(D 



where e M" (respectively, J" e M") contains the 
first (respectively, the second) component of all the C^-'^ 
vectors, for j = l,...,n. Notice that the products by 
matrices A-^, (D'')-^, and (D'")-^ are also implemented 
with O(nlogn) cost via the FFT (instead of the standard 
O(n^) cost of matrix-vector multiplications). For exam- 
ple, A-^v — U^A Uv, where the products by U and U"'^ 
are implemented by the FFT and A is diagonal. 
The proximity operator of .gi(v) = A||v||2, for i = 
1, n, is given by the so-called vector-soft function ^231, 
Il24l . which is the generalization of the soft thresholding 
function for vector arguments: 

prox^llH^(v) — vect-soft(v, a) 

■ max{||v||2 — A, 0}, 



V 2 



with the convention that 0/||0||2 = 0. 
The proximity operator of (7n+i(v) = ^| 
affine function. 



P™X(^/2)l|--y|li(v) 



7y - 



7+1 



V — y||2 IS an 



(15) 



Plugging these equalities into Algorithm |2] we obtain Al- 
gorithm [3] which is very similar to the algorithm presented in 
ll23 1. Of course, the FFT-based implementation in line 6 is only 
vaUd under the assumption that A is a circulant matrix, that is 
that the convolution that it represents uses periodic boundary 
conditions, the same being true for the difference operators 
represented by the matrices. 

Finally, notice that, in general, none of the Dj matrices 
has full column rank, and we also cannot guarantee that A 
has full column rank. However, if A corresponds to a low- 
pass filter with non-zero gain for constant images, matrix 



(Di)^...,(D„)^,A^ 



has full-column rank, and 



the conditions of Theorem 1 are satisfied, thus guaranteeing 
convergence of the ADMM algorithm. 

B. TV-Based Deconvolution with Unknown Boundary Condi- 
tions 

To handle images with unknown boundary conditions, the 
observation model should express the fact that the boundary 
pixels are not observed; that is. 



y = MAx + n, 



(16) 



where n denotes white Gaussian noise of unit variance and (as 
above) A e M" " is a matrix representation of the convolution 
with some kernel; the new element in (fT6] l is M G M™^" 
(with m < n), which is a masking matrix, i.e., a matrix whose 
rows are a subset of the rows of an identity matrix. The role 
of the masking matrix is to observe only the subset of the 



Algorithm 3: ADMM-TV 



1 Set fc = 0, choose a,/3 > 0, Uq^' and Aq \ for 



i = l,...,n + l 
2 repeat 

C ^ Ufe + dfe 

r = [(CW)2,...,(C("))2]^ 



h\2 



10 
11 



12 



13 



/3|A|2)~'u(/?A^C^"+'' +a (D'')^^'' + a (0")^^) 
for j = 1 to n do 



'fe+i 



vect-softlDj Zfc+i — d 



k ' a 



end 



"fc+1 ^ T+:a(^y + PiA Zfc+i - dj. 
flfc+i ^ flfe - (A Zfc+1 - n^^^ ) 



k +~ k + 1 



14 until stopping criterion is satisfied 



image domain in which the pixel values of Ax do not depend 
on the boundary pixels; consequently, the type of boundary 
condition assumed for the convolution represented by A is 
irrelevant, and we may adopt periodic boundary conditions 
for computational convenience. 

Assuming that A models the convolution with a blurring 
filter with a limited support of size {1 + 2 I) x [1 + 2 1), and that 
X € M" and Ax e M" represent square images of dimensions 
^/n X y^, then matrix M e M'"^", with m = - 21)^, 
represents the removal of a band of width / of the outermost 
pixels of the full convolved image Ax. 

Under the model in (fT6l l. TV-based deconvolution with 
unknown boundary conditions is formulated as 



1 

5 = arg min - ||y - M Ax||2 + A^ ||D, 



X 2, 



(17) 



where A and the matrices have exactly the same meaning 
as in dSj. Problem ([TtT i can be seen as hybrid of deconvo- 
lution and inpainting, where the missing pixels constitute the 
unknown boundary. If matrix M is assumed to be the identity, 
problem ( [TtI i reduces to TV-regularized periodic deconvolu- 
tion (as in dHJ). Conversely, ( [TtI i becomes a pure inpainting 
problem if A is identity. Moreover, the formulation ( [TT] ) can 
be applied to problems where not only the boundary, but also 
other pixels, are missing (for example, due to saturation of the 
image sensor). 
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At this point, one could be tempted to map ( fTTT i into (|2]i 
using the same mappings as in (l9])-(fT2b, with M A replacing 
A in ( fT2] i. However, unlike A, matrix MA is not circulant 
(it is not even square), thus it is not diagonalized by the DFT, 
and the inversion in line 4 of Algorithm[3]cannot be computed 
efficiently as in ( fT4b . To sidestep this difficulty, we propose to 
replace the mapping ( fTOl ) by 

5«+i(v) - i||y-Mv||2, (18) 

and leave the mapping ( IT2] | unchanged. The only resulting 
change to Algorithmic] is in line 11, because is no longer 
given as in (fTOb . but is now as defined in ( fTSl l. The proximity 
operator of g„+i(v) = ^||Mv — y||2 is given by 

Pi"0X(-y/2)||M.-y||^(v) = argmin7||Mu-y||^ + ||u- v||^ 

= (7M'^M + I)"i(7M^y + v). (19) 

Notice that ( fTST i is a special case of (fT9] l for M = I. Due to 
the structure of M (its rows are a subset of the rows of an 
identity), matrix M^M is diagonal and so is 7M^M + I, 
its inversion thus being trivial and having 0{n) cost. Observe 
also that M^y is the extension of the observed image to the 
same size as x, by creating a boundary of zeros around it. 

In summary, the algorithm herein proposed for TV-based 
deconvolution with unknown boundary conditions has the 
exact same structure as Algorithm [3] with line 1 1 replaced 
by the following expression: 

u'^lX'^ ^ (M^M + (M^y + /3(A z,+i - 4"+^))) ; 

notice that since (M-^M + /3l) ^ is a diagonal matrix, mul- 
tiplying by it is equivalent to element-wise multiplication by 
a vector. The cost per iteration of the algorithm is dominated 
by the 0{ti log n) of the FFT implementations of the products 
by the DFT matrices U and (line 6). 

IV. Experiments 
A. Deblurring with unknown boundary conditions 

The proposed ADMM approach for image deconvolution 
with unknown boundary conditions was tested on the bench- 
mark images "Lena" and "cameraman", degraded with 4 
different blurring filters (out-of-focus, linear motion, uniform, 
and Gaussian), all of size 9x9 {i.e., (2/ + 1) x (2/ + 1), 
with / = 4), at three different noise levels: 40 dB, 35 dB, and 
30 dB of blurred signal to noise ratio (BSNR). The images 
are blurred using non-cyclic convolution, where the external 
region (of width / = 4) is obtained by copying the boundary 
rows and columns 4 times; that is, the external region above 
the image contains 4 copies of the first row, and similarly for 
the other boundaries. Of course, the deconvolution algorithm 
does not have this information. 

Different ADMM approaches for TV-regularized deconvo- 
lution were considered: (1) ignoring the non-periodic nature of 
the convolution, i.e., assuming (wrongly) that the convolution 
is periodic (this approach will be referred to as "old"); (2) 
using the approach proposed in Subsection IIII-BI of this 
paper (which will be referred to as "new"); (3) preprocessing 



the blurred images using the "edgetaper" Matlab function, 
but ignoring the non-periodic nature of the convolution (this 
approach will be referred to as "tapered"; (4) using the strategy 
described in the appendix, which will be referred to as "CG". 

The regularization parameter A was manually adjusted to 
yield the best ISNR {improvement in signal to noise ratio) for 
each experimental condition. For the ADMM parameters (a 
and P), which affect the convergence speed, we used heuristic 
rules, which lead to a good performance of the algorithm: 
/? was chosen to be proportional to A specifically /3 — lOA; 
concerning a, we used the heuristic rule a = min{l, 5000/3}. 

Tables H] and HI] show the average ISNR values and the 
processing complexities and times of the algorithms consid- 
ered. The processing time is presented in second^ and the 
complexity measured by the number of calls to the FFT 
algorithm required per iteration. The inabiUty of the baseline 
method ("old") to address non-cyclic deconvolution is evident 
from Table H] while the approaches that model the boundary 
as unknown are clearly better. In terms of ISNR, the new 
method proposed in Subsection IIII-BI and the one described 
in the appendix yield essentially identical values, which is 
a natural and reassuring fact, since both methods attack the 
same optimization problem. However, our approach ("new") 
is considerably faster than "CG", having approximately half 
of the complexity per iteration of the method described in 
the appendix. In addition to being faster than "CG", our 
approach is also considerably simpler to implement, since 
it does not involve any inner iterative solver. To serve as 
benchmark. Tables H] and HI] also show results for the baseline 
method ("old") for observed images obtained with periodic 
convolutions. The visual advantage of the proposed approach 
over the "old" and the "tapered" methods are illustrated in 
Figures |2] and |3] The boundary ripple artifacts resulting from 
an incorrect convolution model are quite evident in "old" 
approach and still visible in the "tapered" method. In contrast, 
these boundary artifacts are essentially absent with the "new" 
approach, whose results are visually similar to the (ideal) ones 
obtained for cyclic degradations. The restored images obtained 
with the "CG" method are visually indistinguishable from 
those obtained with the new method, so we don't show them. 

B. Simultaneous deconvolution and inpainting 

The proposed method was also successfully tested in simul- 
taneous non-cyclic deblurring and inpainting. As an example 
of the performance of our method. Figure ]4] shows the results 
obtained for the "cameraman" image, non-periodically blurred 
with a 9 X 9 uniform blur (assuming the same boundary 
conditions explained in the previous experiments), with 40dB 
BSNR, and with 20% of missing pixels. 

V. Conclusions 

We presented a new strategy to extend recent fast algorithms 
for image deconvolution, based on the alternating direction 
method of multipliers (ADMM), to problems with unknown 

^The methods were implemented in Matlab and run on a Intel Core i5 
processor at 2.39GHz. 
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TABLE I 

Comparison of different approaches; the ISNR values are 
computed for the valid pixels and averaged over the two 
images and the four blurring filters considered. 
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TABLE II 

Computational cost of different approaches, measured in 
number of ffts and processing time (seconds) per iteration. 




(a) "old" 



(b) "tapered" 




(c) "new" 



(d) Ideal (cyclic) 



Fig. 2. Deblurring results obtained for tlie "Lena" image blurred with the uni- 
form filter, at 40dB BSNR. From left to right: "old" (ISNR = -2. 67 dBs), 
"tapered" (ISNR = 5.52dBs), "new" (ISNR = 7.09dBs) and "ideal" 
(cyclic degradation, ISNR = 7.52dBs). 



boundary conditions. The proposed approach is also able to ad- 
dress a more general class of problems where the degradation 
model consists of a combination of a convolution (blurring) 
and a loss of pixels. Our approach is validated by adapting a 
recent state-of-the-art algorithm for deconvolution under total 
variation regularization L23J . and comparing it with the the 
recent approach of 11211 . The proposed technique was shown 
to attain results similar to those obtained with the technique 
of II2TII . while being simpler and requiring roughly half of the 
processing time. 

Appendix: TV-based deconvolution with unknown 

BOUNDARY CONDITIONS USING THE REEVES-SOREL 
APPROACH 

A key feature of the approach proposed in Subsection IIII-BI 
is that each step of the ADMM algorithm is computed in 
closed form, without the need for any inner iterations. In this 
appendix, we describe an alternative approach (which does not 
share this desirable property) based on a technique proposed 
in |fT9l , and that was also very recently used in lISTI . The 
goal of this appendix is to adapt the method proposed in 
II2TI to the ADMM formulation and notation used in this 
paper. Let the objective function in ( [TtI i be mapped into the 
form ^ as in (l^-(fT2b, with MA replacing A in ([T2]i. As 
mentioned in Subsection IIII-BI matrix MA is not circulant, 
thus it is not diagonalized by the DFT, and the inversion in 
line 4 of Algorithm [3] cannot be computed efficiently as in 
(fT4l l. Following |21|, that inversion can be addressed using 
the technique introduced in ||T9l . 



With MA replacing A, the matrix to be inverted in each 
iteration of ADMM is given by 

a(D'')^D'' + a(D")^D" + ^A^M^MA; (20) 

although D'\ D", and A are circulant, this matrix cannot be 
inverted as in ( fT4l) . due to the presence of M. As in |fT9l , let 

MA' 
b 

where b contains the rows of A that are missing from MA 
and P is a permutation matrix that puts these missing rows in 
their original positions in A. To simplify the notation, write 
V ^ a(D'')'^D'' + a(D^)^D". Noticing that 



A'^A = [A^M^ , b^jP^P 



MA 
b 



= A'^M^MA + b^b 



(because P is a permutation matrix, thus orthogonal), the 
matrix inversion can be written as 

(V + ;3A^M'^MA) = (V + PA^A - /3h^hy^ 

= C - Cb^(bCb^ - (l//3)I)"^bC, (21) 

where the second equality results from setting C = (V + 
/?A^A)^^ and using the Sherman-Morrison- Woodbury matrix 
inversion identity. The inversion C = (V + /3A^A)^^ ~ 
(ct(D'')'^D'* + ct(D^)'^D^ +/3A'^A)"^ can be efficiently 
computed via FFT, as in (fl4l i. The inner inversion in (|2TI) can- 
not be computed via FFT; however, its dimension corresponds 
to the number of unknown boundary pixels (number of rows 
in b), thus it is usually much smaller than x. As in ||T9l . ||2T1| . 
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(a) •'old" (b) "tapered" 




(c) "new" (d) Ideal (cyclic) 



Fig. 3. Deblurring results obtained for tiie "cameraman" image blurred witii 
the linear motion filter, at 40dB BSNR. From left to right: "old" (ISNR = 
2.03dBs), "tapered" (ISNR = 9.07dBs), "new" {ISNR = W.OSdBs), 
and "ideal" (cyclic degradation, ISNR = 10A7dBs)}. 




(c) Result from (a) (d) Result from (b) 



Fig. 4. Results obtained, with our approach, for the cameraman image 
degraded with a non-cyclic convolution with a 9 X 9 uniform blur, at 40dB 
BSNR, followed by loss of 20% of the pixels: (a) blurred image, (b) blurred 
image with 20% missing pixels, (c) and (d) images recovered from (a) and 
(b), respectively. 



we adopt the conjugate-gradient (CG) algorithm to solve the 
corresponding system; we verified experimentally that taking 
only one CG iteration in each ADMM iteration yields the best 
results and is the fastest option. 
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